# Media coverage and Sentencing 
libraries <- c("lavaan", "mediation", "MASS", "modelsummary", "dplyr", "kableExtra")
for (library_name in libraries) {
  library(library_name, character.only = TRUE)
}
setwd("~/Documents/Journal articles/Identification")  
rm(list = ls())
set.seed(123)
n <- 1000
ep1 <- rnorm(n, mean = 0, sd = 1)
ep2 <- rnorm(n, mean = 0, sd = 1)
ep3 <- rnorm(n, mean = 0, sd = 1)
ep4 <- rnorm(n, mean = 0, sd = 1)
ep5 <- rnorm(n, mean = 0, sd = 1)
X <- exp(ep1) # confounder
D <- pmin(pmax(rpois(n, lambda = exp(0.4*X + ep2)), 1), 13) # media coverage (higher values indicate greater coverage)
linpred_1 <- -0.3*D + 0.4*X + ep3
linpred_2 <- -0.2*D + 0.3*X + ep3 + 0.5  # Added offset to separate predictions
linpred_3 <- -0.1*D + 0.2*X + ep3 + 1.0  # Added larger offset to separate predictions
exp_linpred_m <- exp(cbind(linpred_1, linpred_2, linpred_3))
probs_m <- exp_linpred_m / rowSums(exp_linpred_m)
M <- apply(probs_m, 1, function(p) sample(1:3, size = 1, prob = p)) # public opinion (lesser, neutral, harsher)
linpred_4 <- -0.7*D + 0.4*X + ep4
linpred_5 <- -0.6*D + 0.3*X + ep4 + 0.5  # Added offset to separate predictions
linpred_6 <- -0.5*D + 0.2*X + ep4 + 1.0  # Added larger offset to separate predictions
exp_linpred_b <- exp(cbind(linpred_4, linpred_5, linpred_6))
probs_b <- exp_linpred_b / rowSums(exp_linpred_b)
B <- apply(probs_b, 1, function(p) sample(1:3, size = 1, prob = p)) # judge behavior (lesser, neutral, harsher)
table(M) 
table(B) 
Y <- pmin(pmax(rpois(n, lambda = exp(0.04*D + 0.2 * M + 0.3 * X + 1.8 * B + ep5)), 6), 1500) # sentencing (ranges between 6 to 1500 months)
sem_cim <- data.frame(D, X, M, B, Y)
cma_cim <- data.frame(D, X, M, B, Y)

summary <- cma_cim[, c("D","X", "M", "B", "Y")] %>% 
  select(`Media coverage` = D, `Confounder` = X, `Public opinion` = M, `Judge behavior` = B, `Sentencing` = Y)
summary_list <- lapply (summary, scale)
emptycol = function(x) " "
datasummary(`Media coverage` + `Confounder` + `Public opinion` +`Judge behavior` + `Sentencing` ~  Mean + SD + Min + Max + Heading("Boxplot") * emptycol + Heading("Histogram") * emptycol, data=summary, fmt = "%.3f") %>% column_spec(column = 6, image = spec_boxplot(summary_list)) %>% column_spec(column = 7, image =spec_hist(summary_list))

# CMA (causally independent mediators)
cma_cim$M <- factor(cma_cim$M, levels = c(1, 2, 3), ordered = TRUE)
cma_cim$B <- factor(cma_cim$B, levels = c(1, 2, 3), ordered = TRUE)
model.B <- polr(B ~ D + X, data=cma_cim, method="probit", Hess=TRUE)
model.M <- polr(M ~ D + X, data=cma_cim, method="probit", Hess=TRUE) 
model.Y <- lm(Y ~ D + B + M + X, data=cma_cim)
med.out.B <- mediate(model.B, model.Y, sims=1000, boot = TRUE, treat="D", mediator="B")
med.out.M <- mediate(model.M, model.Y, sims=1000, boot = TRUE, treat="D", mediator="M")
summary(med.out.B)
summary(med.out.M)

                              # SEM Framework
combined_model <- '
  # Regressions
  M ~ c1*D + c2*X
  B ~ c3*D + c5*X
  Y ~ c6*D + c7*M + c8*B + c9*X
  
  # Indirect effects
  IE_M := c1 * c7  # Indirect effect through M
  IE_B := c3 * c8  # Indirect effect through B
  TE := IE_M + IE_B + c6  # Total effect
  DE := c6  # Direct effect
  IE := IE_M + IE_B
'
fit <- sem(combined_model, data = sem_cim, se = "bootstrap", bootstrap = 1000)
sem_cim_summary <- summary(fit, standardized = TRUE, ci = TRUE)
IE_M_estimate <- lavInspect(fit, "est")[["IE_M"]]
IE_B_estimate <- lavInspect(fit, "est")[["IE_B"]]
IE_estimate <- lavInspect(fit, "est")[["IE"]]
DE_estimate <- lavInspect(fit, "est")[["DE"]]
TE_estimate <- lavInspect(fit, "est")[["TE"]]
sem_cim_summary

# CMA (causally dependent mediators)
setwd("~/Documents/Journal articles/Identification")  
rm(list = ls())
set.seed(123)
n <- 1000
ep1 <- rnorm(n, mean = 0, sd = 1)
ep2 <- rnorm(n, mean = 0, sd = 1)
ep3 <- rnorm(n, mean = 0, sd = 1)
ep4 <- rnorm(n, mean = 0, sd = 1)
ep5 <- rnorm(n, mean = 0, sd = 1)
X <- exp(ep1) # confounder
D <- pmin(pmax(rpois(n, lambda = exp(0.4*X + ep2)), 1), 13) # media coverage (higher values indicate greater coverage)
linpred_1 <- -0.3*D + 0.4*X + ep3
linpred_2 <- -0.2*D + 0.3*X + ep3 + 0.5  # Added offset to separate predictions
linpred_3 <- -0.1*D + 0.2*X + ep3 + 1.0  # Added larger offset to separate predictions
exp_linpred_m <- exp(cbind(linpred_1, linpred_2, linpred_3))
probs_m <- exp_linpred_m / rowSums(exp_linpred_m)
M <- apply(probs_m, 1, function(p) sample(1:3, size = 1, prob = p)) # public opinion (lesser, neutral, harsher)
linpred_4 <- -0.7*D + 0.4*X + 0.4*M + ep4
linpred_5 <- -0.6*D + 0.3*X + 0.5*M + ep4 + 0.5  # Added offset to separate predictions
linpred_6 <- -0.5*D + 0.2*X + 0.6*M + ep4 + 1.0  # Added larger offset to separate predictions
exp_linpred_b <- exp(cbind(linpred_4, linpred_5, linpred_6))
probs_b <- exp_linpred_b / rowSums(exp_linpred_b)
B <- apply(probs_b, 1, function(p) sample(1:3, size = 1, prob = p)) # judge behavior (lesser, neutral, harsher)
table(M) 
table(B) 
Y <- pmin(pmax(rpois(n, lambda = exp(0.04*D + 0.2 * M + 0.3 * X + 1.8 * B + ep5)), 6), 1500) # sentencing (ranges between 6 to 1500 months)
cma_cdm <- data.frame(D, X, M, B, Y)

summary <- cma_cdm[, c("D","X", "M", "B", "Y")] %>% 
  select(`Media coverage` = D, `Confounder` = X, `Public opinion` = M, `Judge behavior` = B, `Sentencing` = Y)
summary_list <- lapply (summary, scale)
emptycol = function(x) " "
datasummary(`Media coverage` + `Confounder` + `Public opinion` +`Judge behavior` + `Sentencing` ~  Mean + SD + Min + Max + Heading("Boxplot") * emptycol + Heading("Histogram") * emptycol, data=summary, fmt = "%.3f") %>% column_spec(column = 6, image = spec_boxplot(summary_list)) %>% column_spec(column = 7, image =spec_hist(summary_list))

m.med <- multimed(outcome = "Y", med.main = "B", med.alt = "M", treat = "D", covariates = "X",  data = cma_cdm, sims = 1000)
summary(m.med)

                  # SEM Framework
combined_model <- '
  # Regressions
  M ~ c1*D + c2*X
  B ~ c3*D + c4*M + c5*X
  Y ~ c6*D + c7*M + c8*B + c9*X
  
  # Indirect effects
  IE_M := c1 * c7  # Indirect effect through M
  IE_B := c3 * c8  # Indirect effect through B
  IE_MB := c1 * c4 * c8
  TE := IE_M + IE_B + IE_MB + c6  # Total effect
  DE := c6  # Direct effect
  IE := IE_M + IE_B + IE_MB
'
fit <- sem(combined_model, data = cma_cdm, se = "bootstrap", bootstrap = 1000)
sem_cdm_summary <- summary(fit, standardized = TRUE, ci = TRUE)
IE_M_estimate <- lavInspect(fit, "est")[["IE_M"]]
IE_B_estimate <- lavInspect(fit, "est")[["IE_B"]]
IE_MB_estimate <- lavInspect(fit, "est")[["IE_MB"]]
IE_estimate <- lavInspect(fit, "est")[["IE"]]
DE_estimate <- lavInspect(fit, "est")[["DE"]]
TE_estimate <- lavInspect(fit, "est")[["TE"]]
sem_cdm_summary
